M estimators in R

Nonparametric Statistics
AY 2023-2024

andrea.cappozzo@polimi.it
@AndreaCappozzo
andreacappozzo.rbind.io

Problem setting 1/2

  • Let us try to code the iterative algorithm described in slide 12, but with an additional (robust) scale estimate \(\hat{\sigma}\)

  • Consider the natural logs of the annual incomes of \(10\) people

x <-
  c(9.52, 9.68, 10.16, 9.96, 10.08, 9.99, 10.47, 9.91, 9.92, 15.21)

Problem setting 2/2

plot(x=x,y=rep(0,length(x)), col=c(rep("black",9),"red"),pch=19)
  • Aim: downweight the impact of the red point in \(\hat{\mu}\)

Manual implementation M estimator

manual_M_location <-
  function(x,
           k,
           type = c("Huber", "Tukey"),
           tol = sqrt(.Machine$double.eps),
           itermax = 1000) {
    
    
    mu <- mu_old <- mu_vec <- median(x) # initial value
    rob_scale <- mad(x) # it will be kept fixed
    crit <- TRUE
    iter <- 0
    
    weigth_f <- switch (type,
                        "Huber" = function(x) pmin(1, k/abs(x)),
                        "Tukey" = function(x) ((1-(x/k)^2)^2)*(abs(x)<=k)
    )
    
    while(crit){
      w_i <- weigth_f((x-mu)/rob_scale)
      mu <- weighted.mean(x = x,w = w_i)
      err <- abs(mu-mu_old)
      
      mu_old <- mu
      mu_vec <- c(mu_vec,mu)
      iter <- iter+1
      
      crit <- (err > tol & iter < itermax)
    }
    list(mu=mu, s=rob_scale, w_i=w_i, it=iter)
  }

M estimator with Huber function

fit_M_huber <- manual_M_location(x = x,k = 1.5,type = "Huber")
fit_M_huber$mu
[1] 10.00333
fit_M_huber$w_i
 [1] 0.66717002 0.99731602 1.00000000 1.00000000 1.00000000 1.00000000
 [7] 0.69099748 1.00000000 1.00000000 0.06193319
plot(x=x,y=rep(0,length(x)), pch=19,cex=fit_M_huber$w_i)

M estimator with Tukey function

fit_M_tukey <- manual_M_location(x = x,k = 4.68, type = "Tukey")
fit_M_tukey$mu
[1] 9.959491
fit_M_tukey$w_i
 [1] 0.6547726 0.8516119 0.9221404 0.9999995 0.9715115 0.9981617 0.5513463
 [8] 0.9951664 0.9969210 0.0000000
plot(x=x,y=rep(0,length(x)), pch=19,cex=fit_M_tukey$w_i)